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dhe  correlation  function  of  two  signals  can  rapidly  be  obtained  by 
taking  the  inverse  discrete  fa6t  Fourier  transform  (IDFFT)  of  the  product 
of  the  discrete  fast  Fourier  transform  (DFFT)  of  the  first  signal  times 
the  complex  conjugate  of  the  DFFT  of  the  second  signal.  The  subroutine 
RAPIDC  accomplishes  this  operation. 

In  a preliminary  report  to  the  NUSL,  the  Sperry  Rand  Corporation 
reports  (reference  (a))  success  in  performing  correlation  by  FFT  as 
follows: 

The  discrete  fast  Fourier  transform  (DFFT)  of  signal  1 and  the  DIPT 
of  signal  2 are  taken.  The  complex  conjugate  of  the  latter  is  then 
multiplied  term  by  term,  with  the  former  transform.  The  inverse  discrete 
fast  Fourier  transform  ( IDFFT)  when  properly  scaled  yields  the  correlo- 
gram  of  signal  1 cross  (or  replica)  correlated  with  signal  2. 

Continued  developments  in  the  field  of  the  discrete  fast  Fourier 
transform  (reference  (b))  make  implementation  of  this  concept  easier 
and  more  frequent  (reference  (e)).  The  primary  advantage  of  correlations 
with  FFT  is  the  high  rapidity  with  which  the  computations  can  be  per- 
formed. So  staggering  is  the  time  element  that  it  is  conceivable  to 
disallow  time  domain  replica  and  cross  correlations  by  the  classical  / 
method  for  high  order  cases.  /£ 


Appendix  A is  a program  listing  with  one  applicable  subroutine,  ft  t? 
CO£JAD,  given  in  Appendix  B.  The  subroutine  for*  FFT  can  be  found  in 


ft 
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reference  (b)  or  the  user  may  prefer  to  use  another  transformation  sub- 
routine of  his  own  choosing.  Appendix  C gives  a sample  case  of  cyclic 
and  acyclic  correlation  of  two  pulses. 

THEORY 


Tiie  time  domain  definition  of  the  correlation  function  is  well 
known  for  real  causal  functions,  x(t)  and  y(t  +*f* ) 60  be: 

0«rCr)»^  T \ Kto^t*r)dt 
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Equation  (1A)  from  Bendat  and  Pier sol  (reference  (c))  (digitally 
eaqpressed  as  Eq.  (IB)  in  reference  (i))  is  a classical  definition  re- 
quiring large  and  often  excessive  amounts  of  computation  time  for  large 
arrays  X(N)  and  Y(N)  . The  FORTRAN  implementation  of  this  definition 
for  the  discrete  case  is  straight  forward  and  its  use  throughout  govern- 
ment and  industry  is  quite  common.  The  basic  definition  is  sometimes 
modified  and  a normalized  version  forces  to  fall  between  +1  and 


It  has  long  been  known  that  the  inverse  transform  of  the  cross 
spectral  density  function  is  the  correlation  function.  Davenport  and 
Root  (reference  (d))  state  that: 


jj 
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where  0xy  is  the  correlation  function  and  S^y  is  the  cross  spectral 
density  function.  J 

Bie  cross  spectral  density  function  is,  in  general,  complex  con- 
sisting of  a real  part,  the  co  (coincidental)  spectral  density  function, 
and  an  imaginary  part,  the  quadrature  spectral  density  function  (refer- 
ence (c)). 

(UJ)  of  x(t)  and  y(t) 

e cross  spectrum  is  given 
by: 

Sxy(*^)  = X(cJ)  • Y*(iO)  (3) 

where  indicates  conjugation 

x(t)  X(u)) 

I y(t)  y(oj) 

and  # p indicates  transform  pair. 

Recently,  Singleton  (reference  (b))  has  presented  an  FFT  subroutine 
coded  in  FORTRAN  which  can  be  called  as  follows: 

CALL  FFT(XX,YY,NNN,NNN,NNN,ISN) 

where  (XX  + jYY)  is  a complex  input  and  NNN  is  the  Number  of  Values  to 
be  transformed  (NNN  must  contain  no  prime  factor  greater  than  23)  and 
where  for  ISN  = -1  the  forward  discrete  fast  Fourier  transform  i6  com- 
puted and  returned  as  the  complex  output  (XX  + jYY).  Similarly  for 
ISN  = 1 the  inverse  discrete  fast  Fourier  transform  is  computed  and 
returned  as  the  complex  output  (NNN  times  (XX  + jYY)). 

Given  the  two  real  signals  x(t)  and  y(t)  they  can  be  transformed 
as  the  complex  signal  (x(t)  + jy(t))  yielding  the  complex  transform, 
Z(oi)  , of  the  complex  input  signal.  The  transform  of  x(t)  and  y(t) 
can  be  recovered  as  shown  by  Cooley  (reference  (h))  as  follows:1 


X(W)  = (Z(U»)  + Z*(-W))/2  (If) 

Y(W)  = (Z(W)  - Z*(-W))/2J  (5) 


The  cross  spectral  density  function  S™ 
referred  to  by  Papoulis  (reference  (f))  as  th 


1 A brief  and  complete  derivation  of  Eq. 
Eby  in  reference  (g). 


(4)  and  (5)  is  shown  by 
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Multiplication  term  by  term  of  X*^)  with  Y(to)  can  be  ahown 
from  Eqs.  (4)  and  (5)  to  yield  the  real  (CO)  and  imaginary  (QUAD)  parts 
of  the  cross  spectral  density  function,  S^W)  , as  follows: 

CO(U>)  = 0.50  (XX(OJ)YY(-VaJ)  + XX( - V*))YY(U) ) ) (6) 


JQUAD(U))  = 0.25  (XX(U))2  + YY(U)2  -XX(-W)2  -YY(-Vj))2)  (7) 


where  XX(W)  + JYY(uJ  ) is  the  complex  transform  output  (Z(l*>))  of 
Singleton's  FFT  subroutine. 


Implementation  of  these  two  equations  in  FORTRAN  with  positive  in- 
dexing and  taking  advantage  of  symmetries  present  is  done  by  the  Sub- 
routine COQUAD.  A listing  of  COQUAD  is  enclosed. 


Taking  the  inverse  discrete  Fourier  transform  of  the  CO  and  Quad- 
rature spectral  density  functions,  properly  scaling  the  output  gives 
the  correlation  function  equivalent  to  Eq.  (l). 


Removal  of  the  mean  and  division  by  the  standard  deviation  of  the 
data  prior  to  the  described  procedure  for  RAPID  Correlation  yields  a 
normalized  correlation  coefficient  between  plus  and  minus  one. 

SUBROUTINE  RAPIDC 


Subroutine  RAPIDC  is  called  from  the  main  program  as  follows; 


CALL  RAPIDC ( XX, YY,NUMX,NUMY,N2FL0T, NOWRAP) 


where  XX  is  an  array  (dimensioned  at  least  NUMY  for  NOWRAP  = 1 and 

dimensioned  at  least  NUMX  plus  NUMY  for  NOWRAP  = 0) 

YY  is  an  array  larger  than  or  equal  to  XX  (dimensioned  in 

the  main  program  with  the  rules  above  for  the  XX  array) 

NUMX  is  the  number  of  values  in  XX 

NUMY  is  the  number  of  values  in  YY 


N2PL0T  is  the  number  to  plot  and  is  supplied  by  the  subroutine 

NOWRAP  is  the  flag  parameter  for  cyclic  and  acyclic  correlation; 
NOWRAP  = 0 gives  acyclic  correlation  (i.e.,  zero  wrap 
around)  but  requires  more  array  storage.  See  discussion 
which  follows. 
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CYCLIC  AND  ACYCLIC  OPTION 

In  the  call  to  RAPIDC,  if  NOWRAP  = 0 then  acyclic  correlation  ie 
performed.  However,  in  the  main  program  both  of  the  arrays  must  be 
dimension  to  at  least  NUMX  +NUMY.  For  cross  correlation  of  arrays  with 
8000  values  both  arrays  would  be  dimensioned  to  16000. 

If  NOWRAP  = 1 then  cyclic  correlation  is  performed  and  a wrap 
around  effect  takes  place.  If  a replica  correlation  is  taking  place 
the  user  should  note  that  no  wrap  around  will  occur  until  NUMY  -NUMX 
correlations  t' ke  place.  Hie  value  of  NUMY  -NUMX  will  be  returned  as 
the  number  to  plot  in  cases  where  NOWRAP  =1.  An  example  would  be  of 
the  replica  correlation  of  1000  values  with  16000.  No  wrap  around 
would  take  place  until  15000  correlations  for  NOWRAP  equal  to  1 and 
the  N2PL0T  would  equal  15000. 

The  authors  suggest  as  a rule  of  thumb  that  NOWRAP  be  set  equal  to 
0 for  cross  correlation  with  both  data  arrays  'over'  or  'double' 
dimensioned  and  that  NOWRAP  be  set  equal  to  1 replica  correlations  with 
the  user  only  plotting  N2PL0T  number  of  values  supplied  by  the  sub- 
routine. 

LIMITATIONS  OF  THE  FFT 

Limitations  of  Singleton* 6 FFT  are  that  it  must  be  called  with  a 
highly  composite  number  (that  is  one  with  no  prime  factor  greater  than 
23).  The  subroutine  RAPTDC  calls  FFT  with  NUMY  for  NOWRAP =1  and  calls 
FFT  with  NUMX  + NUMY  for  NOWRAP  =0.  A table  of  highly  composite 
numbers  is  available  in  reference  (b). 

TIMING  CONSIDERATIONS 

Run  and  timed  on  the  UNIVAC  1108  under  EXEC  II  level  22  compiler 
Replica  Correlation  of  1000  values  with  15000  values  took  13  seconds 
with  RAPIDC.  The  exact  same  case  processed  by  the  discrete  version  of 
the  classical  Eq.  (l)  required  182  seconds.  Timing  included  no  input- 
output  manipulations.  This  increase  in  speed  in  the  ratio  of  lk:l  is 
significant  and  should  not  be  overlooked  when  considering  correlation 
computations. 


LOCATION  OF  RAPIDC 

RAPIDC,  FFT  and  COQUAD  are  located  on  File  one  of  CUR  Tape  U2kl. 
Both  the  symbolic  and  relocatable  versions  are  stored  on  U24l,  so  these 
three  subroutines  may  be  used,  compiled,  listed  or  punched  from  tape. 
The  authors  request  that  the  tape  be  brought  into  the  PCF  as  follows: 
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k 


7R  ASG  A = U24l 

57  XQT  CUR 

IN  A 

TRI  A 

The  tape  U24l  should  be  listed  on  a yellow  run  request  card  as  an  input 
tape  (left  column)  and  circled. 

AN  I/O  PACKAGE  FOR  RAPIDC 

An  input-output  package,  computer  program  Sl6l2,  has  been  written 
and  is  stored  on  File  1 of  CUR  Tape  U24l.  Sl6l2  is  a general-purpose 
FFT  replica  correlation  program  with  basically  the  same  input  and  out- 
put as  S1530  (reference  (i)).  The  primary  difference  between  S1530  and 
S1612  is  the  subroutine  which  performs  the  correlation.  S1612  uses 
RAPIDC  while  S1530  used  a classical  time  domain  approach.  Other  coding 
differences  are  required  because  RAPIDC  uses  an  'in  place'  algorithm. 
Still,  run  time  has  been  reduced  tenfold  with  Sl6l2  over  the  previous 
time  domain  method. 
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APPENDIX  A 


SUoKUuUhE  RAPIuC  **  C.R. ARNOLD  ANU  G.C, CARTER  »*  11/OJ/oG 


SUuKOUT  In£  RAPiUC  ( AX  r Y Y • NUMX » NUMY  »NEPLOT  , NO WRAP ) 


COMMENT  - hMPIOL  Is  A RAPID  CORRELATION  SUBROUTINE  bY  FH 


COMMENT  - RAPluc  CkuSS  CORRELATES  X a Wllh  YY 

COMMENT  - mDMY  IS  THE  NUMBER  OF  VALUES  1 .M  ThE  YY  ARRAY 

COMMENT  - nUmX  IS  hit.  NUMBER  OF  VALUES  IN  ThE  XX  ARRAY 

comment  - ihE  xa  array  is  shorter  ok  Equal  ro  the  ty  akkay 
COMMENT  - NUMX  AND  nUMa+NUMY  MUST  at  HIuhLY  composite 
iThAT  IS  — NO  PRIME  FACTOR  GReATlR  THAN  u.i> 

COMMENT  - i.EPLOT  Is  the  NUMoeR  of  values  to  BE  PLOTTED  bY  The  user 

comment  - NfcPLOl  IS  SUPPLIED  ro  THE  USER  BY  THE  SUBROUTINE 

COMMENT  - nDwRAP  Is  A FLAG  PARAMETER  FOR  CYCLIC  AND  ACYCLIC  CORRELATION 

(THAT  IS  — THERE  IS  NO  WRAP  AROUnU  IF  NOlnRAP=G  — HOWEVER 
Mure  sTuRAGE  IS  kEQUIRED  TO  SUPPRESS  wrap  around) 

COMMENT  - aF  TER  PRuCeSSInG  aX(I)  HOLDS  The  COKRELOGRaM 
COMMENT  - rOR  N0«RaP-1»  BOTm  XX  ANU  YT  MUST  bE  DIMENSIONED  To 
i.UMY  In  IhE  MAIN  PROGRAM 

COmMEnT  - F uk  NUwKaP-0 t UOlH  XX  ANu  YY  MUST  BE  DIMENSIONED  TO 
nUmX  + NuMY  IN  The  main  PROGRAM 
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SPECIFICATION  ANu  TYPt  STATEMENTS 


ulMtNSION  AXU)'YY(l) 


TEST  TO  INSUkE  NumY  IS  GREATER  THAN  OR  EQUAL  TO  NupX 

IFinuMY.GE.NcMA)  Gc  TO  200 
PRINT  ISO 

iso  format i ini i //ioa* incorrect  arguments  three  and  four  in  raPIDc*//) 
3T0P  RakIlic 


CALCULATE  CO.xSTAi.TS 

200  nUMXPI=NUMa+1 
NUMYPI=NUMi+i 
NUMXPY=NUMaaNUMY 


TEST  HO*  TO  i- ILL  arrays  *ITH  ZEROS 
IF inorKAP) juO» ASO>OoO 


IF  Nc»RAP=l  CALCULATE  WITH  WRAP  AKOUNo 

JOu  n2PcOT=nOM(-NUMa 

jO  A00  ISNuMXPl  r NUtsY 
xX(i)=o.O 
AOO  CONTINUE 

nE*NnN=NuM 1 
uu  TO  600 


IF  Nu»RAP=U  CALCULATE  rITH  ZtKO  *RAP  AROUND 

comment  - AA  A No  Y T arrays  MUST  tJE  uIMEnsIONEO  NUMAPY  - i.OTE  - 
4S0  n2PL0T=N0Mi 

uO  3JU  1=NumaP1«NUrXPY 
aXU)-U,0 
30  U cOnIINUE 

oO  330  1=Nj«YP1.n0mxPY 
YYU)=0.O 
SSo  CONTINUE 

nErNNNenOMaP 1 
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CGMPjTt  INVERSE  FAST  FOURIER  TRANSFORM  BY  SlNSLETON'S  METhOU 

1Si«=1 

(.All  FFl  lXA»YY*NE*NiMN»Ntl*>'thi't»NE*NNNfISN) 


SCALu  INVcKSc.  FFT  OUTPUT 

p LIVUM  a— FlO  a 1 (NUwA) 

FLIVE«N=FLOaI  INEcNNN) 
FSCALE=FLNoFiX*FlNE*N 

IF  (NOwRAP)  /1u»  /OOr  no 
710  -0  720  1=I«NUMY 

aX(I)=XaUJ/FSCaLE 
720  lOnTINUE 
oO  TO  7S0 

700  uU  710  1=1 »nuMXPY 
aX(1)=XXU>/FSCa(.E 
710  LOnIINUE 
7S0  oONUNUt 


c 

c Kt TURN  TO  MAIN  PROGRAM  11TH  CORRELOSKAM  IN  XX 
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«t TURN 


C.NO 
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uiMtNSlON  AAUJrYYIl) 
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C 

1HHP1=  ( Nfcriiii«N/ii ) +1 
C 

c 

AXU>=XXU)*YYU) 

YY11)=0.U 

c 

uO  67b  K=*rlhiNPi 
j=rtNPi£”K 

AXX=U.bG*UAlK>*YY< j)+AAlJ).YY(K)  ) 
YY(K)=U.«!5*lAX(J)**i;+YY(  J)  »*2-XX (K ) **2-Y  YIK)**2) 
YY(J)=-YY(a) 

AX(X)=XXX 

XX(J)=XA(K) 

67b  (.unTlNUt 
C 
C 
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